///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/diis.h
///
/// \brief   header file for DIIS acceleration
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_DIIS_HEADER
#define FORMIC_DIIS_HEADER

#include<vector>
#include<algorithm>

#include<boost/shared_ptr.hpp>

#include<formic/vector/vector.h>
#include<formic/vector/matrix.h>

namespace formic {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  uses diis to produce an improved solution based on a history of parameter and error
  ///         vectors
  ///
  /// \param[in]      pmt_hist      the history of the parameters
  /// \param[in]      err_hist      the history of the error vectors
  ///
  /// \return  the improved parameter vector
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  formic::Vector<double,int> diis_accelerate(const std::vector<formic::Vector<double,int> > & pmt_hist,
                                             const std::vector<formic::Vector<double,int> > & err_hist) {

    //// print a greeting
    //if (formic::mpi::rank() == 0)
    //  formic::of << "performing DIIS acceleration." << std::endl << std::endl;

    // ensure history vectors are sane
    if ( pmt_hist.size() != err_hist.size() )
      throw formic::Exception("history vectors are different sizes in formic::diis_accelerate");
    if ( pmt_hist.size() < 2 )
      throw formic::Exception("history length must be at least 2 in formic::diis_accelerate (it is currently %i)") % pmt_hist.size();

    // compute the B matrix
    const int bdim = pmt_hist.size() + 1;
    std::vector<double> bmat(bdim*bdim, 0.0);
    for (int i = 1; i < bdim; i++) {
      bmat.at(0+bdim*i) = -1.0;
      bmat.at(i+bdim*0) = -1.0;
    }
    for (int j = 0; j < bdim-1; j++)
    for (int i = 0; i < bdim-1; i++)
      bmat.at((i+1)+bdim*(j+1)) = formic::xdotc(err_hist.at(j).size(), &err_hist.at(i).at(0), 1, &err_hist.at(j).at(0), 1);

    // prepare the vector for the linear equation
    std::vector<double> bvec(bdim, 0.0);
    bvec.at(0) = -1.0;

    // compute the coefficients for the linear combinations of past parameter vectors
    {
      int info = 0;
      std::vector<int> ipiv(bdim, 0);
      formic::xgesv(bdim, 1, &bmat.at(0), bdim, &ipiv.at(0), &bvec.at(0), bdim, info);
      if ( info != 0 )
        throw formic::Exception("formic::xgesv failed with info = %i in formic::diis_accelerate") % info;
    }

    // compute the improved parameter vector
    formic::Vector<double,int> imp_pmt(err_hist.at(0).size(), 0.0);
    for (int i = 1; i < bdim; i++)
      imp_pmt.axpy(bvec.at(i), pmt_hist.at(i-1));
    //imp_pmt.axpy(-1.0, pmt_hist.at(0));

    // return the improved parameter vector
    return imp_pmt;

  }

}

#endif
